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Abstract 

Critical properties of lattice gases with nearest-neighbor exclusion are investigated via the 
adaptive- window Wang-Landau algorithm on the square and simple cubic lattices, for which the 
model is known to exhibit an Ising-like phase transition. We study the particle density, order 
parameter, compressibility, Binder cumulant and susceptibility. Our results show that it is pos- 
sible to estimate critical exponents using Wang-Landau sampling with adaptive windows. Finite- 
size-scaling analysis leads to results in fair agreement with exact values (in two dimensions) and 
numerical estimates (in three dimensions). 
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I. INTRODUCTION 



Recently efficient methods for estimating the number of configurations of classical sta- 
tistical models have been developed. If the number Q(E) of configurations with energy E 
is determined to sufficient accuracy, many thermodynamic quantities can be obtained with 
little further effort, for any desired temperature. Wang-Landau sampling (WLS) promises to 
be a simple and reliable approach for estimating £l(E) using Monte Carlo simulations 5, 6]. 
Critical exponents, the transition temperature (or chemical potential) and other quantities, 
such a cumulants, may then be estimated via finite size scaling (FSS) analysis lH4(. In this 
paper we test the Wang-Landau algorithm with adaptive windows applying it to the 
lattice gas with nearest-neighbor exclusion. 

Lattice gases have been used extensively as models of simple fluids, and along with the 
Ising model have received much attention in equilibrium statistical physics as a prototype 
for phase transitions. A particularly simple case is the lattice gas with nearest-neighbor 
exclusion (NNE), corresponding to an interparticle potential that is infinite for distances 
< 1 (in units of the lattice constant) and zero otherwise. In the absence of an energy scale, 
temperature is not a relevant parameter, and the system is termed athermal. It is known that 
on bipartite lattices, the lattice gas with NNE suffers a continuous phase transition between 
a disordered phase and an ordered one at a critical value of the density or of the reduced 
chemical potential fi = (311 [1CH13|. (p denotes the chemical potential.) In the ordered phase 
the occupation fractions of the two sublattices are unequal. The grand partition function is 



N„ 



E(z,L)= ^2z»n(N,L), 



(1) 



N=0 



where z = e M is the fugacity, N max is the maximum possible number of particles, and Q(N, L) 
the number of distinct configurations with iV particles satisfying the NNE condition, under 
periodic boundaries. (On a hypercubic lattice of L d sites in d dimensions, N max = L d /2 for 
L even.) The order parameter is the difference between the occupations of sublattices A 
and B: 



xeA xeB 

where <r x is the indicator variable for occupation of site x. 



(2) 



The NNE lattice gas has been studied on various structures: the square llMl3| , triangular 
simple cubic 10], hexagonal 15 L body-centered cubic 16 1, and face-centered cubic 



lattices 



17| . and in higher dimensions 18|. Repulsive lattice gases with exclusion extending 



to second or further neighbors have also been studied |13l . |19I |. Various techniques have 
been applied to study the phase transition, including exact enumeration, series expansion 
(high- and low-density expansions), the cluster- variation method, transfer matrix analysis, 
and Monte Carlo simulation. The model exhibits Ising-like universality on the square and 
honeycomb lattices, while on the triangular lattice (Baxter's hard-hexagon model 20j, |2l|) 



it belongs to the three-state Potts model universality class. 

In this paper we calculate the critical properties of the NNE lattice gas on the square 
and simple cubic lattices using the adaptive-window Wang-Landau (AWWL) algorithm, 
which has been shown to improve the performance of WLS 0|. The critical density 
and chemical potential, as well critical exponents and the reduced fourth-order cumulant, 
are estimated using FSS analysis. The balance of this paper is organized as follows. In 
Sec. II, the adaptive-window Wang-Landau algorithm is briefly reviewed. Sec. Ill contains 
our results for the number of configurations (exact enumeration and simulation results), 
thermodynamic quantities and critical exponents. A summary is provided in Sec. IV. 



II. METHOD 



Consider a statistical model with a discrete configuration space, and let $ denote a 
variable (or set of variables) characterizing each configuration, such as energy or particle 
number. For a given system size, knowledge of the number of configurations (called the 
"density of states" ) for all allowed values of $ permits one to evaluate the partition function 
and associated thermal averages for arbitrary values of the temperature. Wang-Landau 
sampling (WLS) furnishes estimates of the configuration numbers, which we denote 

by reserving to denote the exact values, which are in general unknown. This is 

done by performing a random walk in configuration space, with an acceptance probability 
proportional to 1 / 0(^9') , where denotes the values associated with the newly generated 
(or trial) configuration. In WLS one aims for equal numbers of visits to each set of allowed 
values of as reflected in the histogram, H(i9). 

Various strategies have been proposed to improve WLS and optimize its convergence 



BBS 



-25]. In this work we apply one such scheme, adaptive- windows WLS (AWWLS). This 
method estimates the density of states by determining the range over which the histogram 
has attained the desired degree of uniformity at various stages of the simulation. 

The AWWLS procedure furnishes estimates, Q(N\L), of the number of A-particle con- 
figurations on a lattice of L d sites, to within an overall multiplicative factor which is in- 
dependent of N. For each accepted A-particle configuration, we update the histogram: 
H(N) — > H(N) + 1. Since the density of states is not known a priori, we set Q(N) = 1 
for all N, at the beginning of each sampling level. In the simulation, if N and N + AN 
are the particle numbers in the current and trial configurations, respectively, (in practice, 



AN = ±1), then the acceptance probability is 

tt(N) 



p(N -> N + AN) = min 



n(N + AN) 



1 



(3) 



(To simplify the notation we suppress the dependence of Q on system size L.) 

Whenever a move to a configuration with N particles is accepted, the density of states 
Q(N) is updated, multiplying it by a modification factor / > 1, so: Q(N) —> f ■ Q(N). If 
the trial configuration is rejected we update Q(N) (as well as the histogram) of the current 
N value. As is usual in WLS, the modification factor is initially set to fo = e = 2.71828.... 
After m = 10 4 Monte Carlo steps we check if the histogram satisfies the flatness criterion 
on the minimal window, of width W = (N max — N m i n )/n, beginning with N m i n j^]. The 
histogram is said to be flat if, for all levels in the window of interest, H(N) > 0.8H, where 
the overline denotes an average over levels within the proposed window. If it is not flat, we 
perform an additional m Monte Carlo steps and check again, repeating until the histogram is 
flat on the minimal window. Once this condition is satisfied, we check whether the histogram 
is flat on a larger interval. Thus we define one window and repeat the procedure on the rest 
of the range of N values, forming windows for each stage of sampling. The window positions 
depend on the portion of the histogram that is flat; we include an overlap of three levels 
between adjacent windows. This process is repeated until all values of N have been included 
in a window with a flat histogram. Then a new stage is initiated: the modification factor 
is reduced, / — > v7; an d we reset H(N) = for all N. This process is iterated and the 
simulation halted when / — 1 is approximately 10~ 7 . As explained in 7J, window boundaries 
are not allowed to take the same positions on subsequent stages, to avoid distortions in Q 



that arise when using fixed windows. 

For simplicity, two kinds of trial moves are employed: insertion and removal of particles. 
Including particle-displacement moves - for which N does not change - leads to the same 
results to within uncertainty. We use the R1279 shift register random number generator jg]. 



III. RESULTS 

We study the hard core lattice gas defined above using AWWLS. The estimates Q(N) 
are used to calculate (N) and var(iV) directly; other thermal averages are given by 

Here (A) N is the microcanonical average of quantity A over all configurations having ex- 
actly iV particles, which must also be estimated during the simulation. The development of 



reliable methods for estimating microcanonical averages is an important open problem [28]. 
In the WL procedure, all configurations having the same N should occur with the same 
probability, so that, in principle, the microcanonical average (A) N should be taken over all 
accepted configurations having exactly N particles, with equal weights. We nevertheless 
obtain better results if we restrict the microcanonical averages to the later stages of the 
sampling. Specifically, the averages (4>) N , (<p 2 ) N and (cj) A ) N calculated using all A-particle 
configurations accepted during the simulation yield estimates for critical exponents that de- 
viate significantly from their expected Ising model values. Such deviations are not observed 
for systems with L < 100, but do appear for larger sizes. Similar problems were found in 



studies of spin models using WLS [42|. The results for microcanonical averages improve 
when we restrict the sample to configurations accepted in the later stages of the simulation, 
i.e., for / < 1 + 10- 4 . 



A. Transfer-matrix analysis 

As a preliminary test of our method, we compare our simulation estimates, fl(N), with 
the results of an exact enumeration of Q(N), on a lattice of 8 x 8 sites. The latter are 
obtained via a transfer matrix approach. One begins by enumerating the allowed con- 
figurations {ci, Cm} on a ring of L sites, and storing the number n(cj) of particles in 



each configuration. An A4 x Ai matrix T is then constructed, with T(c J ,Cj) = 1 if adja- 
cent rings may assume configurations q and Cj without violating the NNE condition, and 
T(q, Cj) — if the condition is violated. Then the allowed configurations on an L x L lattice 
with periodic boundaries are those sequences {ci, c 2 , ci} of ring configurations satisfying 
T(ci, C2)T(c2, C3) ■ • • T(cl_i, Cl)T(cl, Ci) = 1. The resulting numbers of configurations for 
L = 8 are listed in Table [H 

To compare our simulation estimates against the exact enumeration, the former must be 
normalized, as simulation in fact provides V(N) = aQ(N), with a an unknown constant, 
independent of N. To eliminate a we multiply the simulation estimates by a factor A, 
varying A so as to minimize ^Ar[Ar(iV) — fl(N)} 2 . This procedure is applied to each of the 
fifteen independent simulation studies, leading to the estimates and uncertainties listed in 
the Table. The relative error in estimating lnf2 is very small except near the minimum and 
maximum occupations. Even in the worst case, N = N max = 32, the relative error in 
is ~ 0.6%. If sampling errors were restricted to these regimes for larger system sizes, the 
effect on estimates for critical properties would be negligible. 

B. Square lattice 

In this work we study 18 system sizes in the range 16 < L < 256 using AWWLS. 
Let N*(L) be the value of N that maximizes £l(N), and let N C (L) be the value of N 
that maximizes the probability distribution P(N) = Q(N) exp\p c N] at the critical point 
H c . Preliminary studies on lattices with L < 200 reveal that N* ~ 0.227L 2 , whereas 
N c ~ 0.369L 2 , so that N c 3> N*. As a result, configurations with N < N* make a negligible 
contribution to thermal averages in the neighborhood of the critical point. To economize 
processor time we therefore restrict our high-statistics studies (which extend to L = 256), 
to iV values between N*(L) and N max = L 2 /2. For L = 8, a study with sampling restricted 
to N > N* = 14 yielded results of equal accuracy as those obtained using unrestricted 
sampling, when compared against exact enumeration. 

The following observation suggests that a further economy of processor time could be 
realized in studies of the critical region. Let P C (N, L) = z?Q(N, L) be the contribution 
to the grand canonical partition function due to the set of all iV-particle configurations 
at the critical point, and let P*(L) = maxAr[P c (iV, L)\ = P(N C ). For N values such that 
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TABLE I: Comparison with numerical results (L=8) for the density of states. The relative error 
is e{N) = (lnO(iV) - ln^(iV))/lnO(vY). 

P C (N , L) / P* (L) < 10~ 50 , say, the contribution to H and thermal averages is negligible. 
It therefore seems reasonable to restrict the sampling to the interval [Ni,N 2 ] of iV values 
such that P C (N,L)/P*(L) > lO" 50 . In practice we use iVi = [0.315L 2 - 131] and N 2 = 



[0.423L 2 + 129], where the brackets denote the largest integer. For L < 300, the largest 
size considered, the resulting interval is small enough to be studied using WLS without 
windows. Surprisingly, restricting the sampling in this manner yields estimates for critical 
exponents that deviate significantly from their expected Ising model values (for example, we 
find j/v — 2.02(5) on the square lattice). We conclude that restricted sampling distorts the 
estimates for the numbers of configurations, and adversely affects microcanonical averages. 

We turn now to the results obtained using the sampling interval [N*(L), N max (L)}. Figure 
[1] shows the number of configurations versus density p; the very good data collapse confirms 
the expected scaling 

n(N,L)~exp[L d g(p)] (5) 

The inset shows N* as a function of system size. Here and below all averages and uncer- 
tainties are obtained using fifteen independent runs. 
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FIG. 1: (Color online) Square lattice: lnfi/L 2 versus density, system sizes as indicated. Inset: 
N*(L) versus system size. 

Two thermodynamic properties used to characterize the transition in lattice gases are the 



particle density p(/x) and the compressibility, 



l 2 ((p 2 ),-(p);) 

(P)l 



(6) 



Figure [2] shows simulation results for the density p and compressibility k as functions of the 
chemical potential; Fig. [3] shows the order parameter, Eq. fl2]) and the susceptibility, 



xM = £ 2 «0 2 ) - (4)1). 



(7) 



The insets in Fig. 3 show the data collapse obtained using the exact critical exponents, 
7/V = 7/4 and (3/v — 1/8, and the high-precision result for the critical chemical potential 



obtained by Guo and Blote 
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ll|, p c = 1.33401510027774(1). 
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FIG. 2: (Color online) Square lattice: density versus chemical potential on the square lattice. 
Error bars are smaller than the symbols. The inflection close to the transition point is weak, being 
imperceptible for the smaller systems. Inset: Compressibility versus chemical potential. 



yzed the dimensionless ratio Q± = (4> 2 ) 2 / (0 4 ) (Fig.HJ), related to Binder's reduced 



We ana 

cumulant 29j, which is expected to take a universal value at the critical point. Let fi c ,Li 
denote the chemical potential at which Q^Li) = Qi(L i+1 ) = Q c ,l^ i-e., the crossing between 
cumulants associated with a pair of successive systems sizes Li and Li + i, and let L = 



0.8 



0.6 



0.4 



0.2 








1500 



1200 



900 



600 



300 







L=48 
L=72 
L=128 
L=196 
L=256 




I 

I 

if 



-100-50 50 100 



I 

it 

In 



I Si 



<XKXXXXXXX), 



0.9 



1 



1.1 1.2 



1.3 



1.4 



FIG. 3: (Color online) Square lattice, upper panel: order parameter versus chemical potential. 
Lower panel: susceptibility versus chemical potential. The insets are data-collapse plots. 



y/LiLi+i denote the geometric mean of two successive sizes. It is common to plot [i c ^ and 
Q c ,l versus l/l} v to estimate the critical chemical potential and cumulant, via extrapolation 
to L — > oo. In the present case, however, we observe no tendency; all values of Q c ,l and 
fi c ,L agree to within uncertainty. Averaging over all values we obtain /i c = 1.335(3) and 



Q c = 0.852(6). Though of low precision, these results are consistent with the literature 
values quoted in Table [TTJ 
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FIG. 4: (Color online) Square lattice: fourth-order cumulant for system sizes as indicated. 



Estimates for the critical chemical potential \i c are obtained via analysis of the chemical 
potential values associated with the maxima of the susceptibility and compressibility for 
each system size. The extrapolated values, p c = 1.330(1) using the susceptibility and fi c = 
1.337(2) using the compressibility, are obtained using the susceptibility data for L = 20 - 196 
and the compressibility data for L = 26 - 196. Pooling our results, we obtain p c = 1.332(2). 
Linear extrapolation of the density p(p c , L) versus 1/L yields p c = 0.3675(5) (see Fig. [5]), 
inset), consistent with the critical density reported in [ll|, p c = 0.3677429990410(3). (Using 
the precise estimate for p c quoted above 11] . we obtain p c = 0.36800(5).) 

Applying FSS analysis to the results for susceptibility for L = 22 - 256 yields j/u — 
1.764(7) (see the inset of Fig. [9]). To estimate j3/u we analyze <f) c (L) using the above cited 
value of p c [ill ]: our data for L < 256 yield f3/v = 0.123(2); (see Fig. [TUl inset). (Using 
our own less accurate estimate, p c = 1.332(2), we obtain f3/v = 0.130(9).) Figure [5] shows 
the maximum of the compressibility versus system size; the results are consistent with 
ftm,L ~ InL, as expected for a model in the 2d-Ising universality class. Table HT1 summarizes 
our main results. It is interesting to note that we obtain essentially the same results, 
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FIG. 5: (Color online) Square lattice: maximum of the compressibility versus InL. The inset shows 
the density ELS £1 function of 1/L ' v . 

7/V = 1.764(8), and (3/v = 0.122(3), if we exclude the data for the two largest system sizes 
from the analysis. 
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TABLE II: Critical values for the square lattice obtained via WLS with adaptive windows. The 
results from [31[ were obtained using Monte Carlo simulations while Refs. ll|, |3Q . |32j use a 
transfer- matrix technique. 



C. Cubic lattice 



We apply AWWLS to the NNE lattice gas on the simple cubic lattice, in system sizes 
L = 8, 10, 12, 14, 16, 18, 20, 22, 24, 28, 32, 40 and 48. In this case we sample the full range of 
iV values, using adaptive windows as described above; for L = 48, we use approximately fifty 
windows. Figure [6] shows Q(N), again verifying Eq. (JSJ). The density and compressibility 
are plotted versus chemical potential in Fig. [TJ while Fig. E] shows the order parameter and 
susceptibility, and the inset of Fig. |6]the fourth-order cumulant. 




FIG. 6: (Color online) Simple cubic lattice: O(iV) versus density, system sizes as indicated. Inset: 
Fourth-order cumulant. 



Using the critical exponent v = 0.6301(4) 36J we plot the values p c (L) (corresponding to 



the maxima of the susceptibility and the compressibility) versus l/L 1 ^. Extrapolation of 
the data for L > 14 yields \x c = 0.05516(9) using the susceptibility, while the compressibility 
data (for L > 24) yield \x c = 0.0567(2). (It is not surprising that the critical value obtained 
using the compressibility is less precise than that found using the susceptibility, as the former 
exhibits a weaker singularity than the latter.) Using our best estimate fi c = 0.05516(9) we 
calculate p(/x c ,L); linear extrapolation (for L > 24) versus l/L 1 /" yields p c = 0.21082(5) (If 



we instead use the estimate fi c = 0.05443(7) [10|j, we find p c = 0.21058(5)). 



Proceeding as in the case of the square lattice, we estimate the critical chemical potential 
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FIG. 7: (Color online) Simple cubic lattice: density versus chemical potential. Inset: Compress- 
ibility versus chemical potential. 
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FIG. 8: (Color online) Simple cubic lattice: Order parameter versus chemical potential. Inset: 
Susceptibility versus chemical potential. 

fi c and the critical moment ratio Q c . Plotting fi cL against we obtain fi c = 0.0552(7) 



via linear extrapolation. This is consistent with previous results lOj which found fi c = 
0.05443(7). Extrapolation of Q c ,l as a function of yields Q c = 0.652(5), somewhat 

higher than the literature value (see Table [TTT]) . If we use \i c = 0.05443(7) 10] to calculate 
Ql{^c)i we observe no significant dependence on L; averaging over all values for L > 18 
yields Q c = 0.636(3). 



FSS ana. 
0.05443(7) 



ysis of the susceptibility furnishes j/v — 2.056(6) (Fig. 15]). Using \i c 



10j . FSS analysis of the order parameter at fi c yields (3/u = 0.504(8) (Fig. 



TDD. (Using our own best estimate, fi c = 0.05516(9), we find (3/v = 0.477(7)). Table HI] 
summarizes our principal results for the simple cubic lattice. As in the case of the square 
lattice, our results do not change significantly if we exclude the data for the two largest 
system sizes from the analysis. 




FIG. 9: (Color online) Maximum of susceptibility versus system size on the square (inset) and 
simple cubic lattices. The solid lines are linear fits used to estimate 



D. Critical behavior of g(p) 



As is known 



lOj, the compressibility of the NNE lattice gas diverges as k ~ a in the 



vicinity of the critical point. (Here /i = (/i — fi c )/^c is the reduced chemical potential.) On 
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FIG. 10: (Color online) Critical order parameter versus system size on square (inset) and simple 
cubic lattice. The solid lines are linear fits used to obtain f3/v. 
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TABLE III: Critical values for the simple cubic lattice obtained via adaptive-window WLS. The 
results from 10(] were obtained using a cluster algorithm of the NNE lattice gas. Refs. 33J and 
34| are from high resolution cluster simulations of the Ising model. 



the other hand it is easy to show that k oc l/\g"(p)\, where g" denotes the second derivative 
of g (denned in Eq. (jSJ) with respect to p. Thus the singularity in the compressibility is 
reflected in a singularity in g. While plots of g versus p appear quite smooth (see Figs. [Hand 
[6]), the second derivative does indeed exhibit a singularity near p c . To obtain g" we perform 
quadratic fits to g(N) on windows of b successive iV values. We choose b large enough 



to eliminate small-scale fluctuations, but small enough that the singular behavior remains 
evident {40]. Despite the rounding incurred by such averaging (in addition, of course, to 
finite-size rounding), the data shown in Fig. [11] provide a clear indication of a developing 
singularity. The figure also shows that the minimum of \g"\ appears to approach zero as 
L — > 00. Our data, however, are not sufficiently precise to verify the expected scaling, 

£,-a/". 
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FIG. 11: Second derivative of g(p) on the simple cubic lattice, system sizes L = 12, 16, 24, 32 and 
48. The inset is a similar plot for the square lattice for system sizes L = 24, 48, 128 and 196. 



IV. CONCLUSIONS 



We perform adaptive-window Wang-Landau simulations of the lattice gas with nearest- 
neighbor exclusion on the square and simple cubic lattices. On the square lattice, comparison 
with an exact enumeration of configurations for L = 8 yields excellent agreement. Using 
finite-size scaling analysis of data for systems of up 256 2 sites on the square lattice and 48 3 
sites on the cubic lattice, we estimate the critical exponent ratios 7/V and /3/v, and the 
critical values of the chemical potential, the density, and the fourth cumulant. In general, 
fair agreement is observed with literature values. In three dimensions, the critical point is 



obtained with an error of about 1.3% compared with previous studies, and exponent ratios 
and Q c with an error of about 3.5%. The precision of our results is considerably less than 
that obtained using transfer-matrix or high-resolution Monte Carlo simulations. Although 
this is somewhat disappointing, we note that despite the widespread interest in Wang- 
Landau sampling, few studies have been published in which critical exponents are obtained 



via this technique 
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LN45| . Our results confirm that it is possible to obtain reasonably 
accurate values for critical exponents and related quantities using Wang-Landau sampling 
with adaptive windows. It thus appears worthwhile to seek further improvements in the 
method, in efforts to develop a simple and versatile approach for studying phase transitions 
via Monte Carlo simulation. 

In this regard two observations seem pertinent. The first is that restricting sampling to 
a subset of densities may worsen the results, even though densities outside this subset make 
a negligible contribution to thermal averages. It appears that the imposition of reflecting 
barriers on the random walk in configuration space somehow distorts the sampling. The 
second point is that the quality of the results furnished by the WLS procedure appears to 
decay with increasing system size, even while maintaining the same flatness criterion and 
schedule of updates of the factor /. Thus, including larger systems sizes in the analysis may 
not improve results; more extensive sampling of small and intermediate system sizes may 
represent a more effective allocation of computing resources. We hope to explore the reasons 
for, and implications of these observations in future work. 
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